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Abstract 

This paper studies growth, percolation, and correlations in disordered fiber 
networks. We start by introducing a 2D continuum deposition model with ef- 
fective fiber-fiber interactions represented by a parameter p which controls the 
degree of clustering. For p = 1, the deposited network is uniformly random, 
while for p = only a single connected cluster can grow. For p = 0, we first 
derive the growth law for the average size of the cluster as well as a formula 
for its mass density profile. For p > 0, we carry out extensive simulations on 
fibers, and also needles and disks to study the dependence of the percolation 
threshold on p. We also derive a mean-field theory for the threshold near 
p = and p = 1 and find good qualitative agreement with the simulations. 
The fiber networks produced by the model display nontrivial density correla- 
tions for p < 1. We study these by deriving an approximate expression for 
the pair distribution function of the model that reduces to the exactly known 
case of a uniformly random network. We also show that the two-point mass 
density correlation function of the model has a nontrivial form, and discuss 
our results in view of recent experimental data on mass density correlations 
in paper sheets. 

Key words: Continuum percolation, deposition models, fibre networks, spatial correla- 
tions 
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I. INTRODUCTION 



There are many phenomena in nature that can be viewed as deposition processes where 
various transport mechanisms bring particles to a surface. These include many different 
examples such as deposition of colloidal, polymer and fiber particles Many deposition 

phenomena involve particles whose size is large compared to the their mutual interaction 
range, and so the main deposition mechanism is due to particle exclusion. Among the 
most studied in this class are the Random and Cooperative Sequential Adsorption models 
There particles are deposited on a surface and either stick or are rejected according 
to certain exclusion rules, with a maximum coverage (the "jamming limit") less than unity. 
These types of models should be contrasted to the case of multilayer surface growth [|l]||[7|], 
where the main focus is in the asymptotic behavior of the growing surface in the continuum 
limit f§. 

A particularly interesting and challenging example involving particle deposition can occur 
in the case of colloidal suspensions. For some such systems, the interparticle repulsion is 
strong enough to prevent multilayer growth [[|. However, the existence of dispersion forces 
can cause the particles to flocculate, or aggregate, and to precipitate out of the suspension 
0,0]. For larger particles or clusters of particles, gravity often induces sedimentation out of 
the suspension. Experiments reveal that sedimentation produces nontrivial spatial structures 
IIJ. A full microscopic treatment of sedimentation is a formidable task, however ]TT|. 

An interesting practical deposition process is the formation of laboratory paper. Paper 
is manufactured from a colloidal suspension, out of which fibers filtrate onto a wire mesh, 
leaving behind a fiber network [[|. Recent experimental measurements of mass density 
correlations in laboratory made paper sheets reveal |12| that the spatial distribution of 
fibers may not be uniformly random. In particular, for low mass density paper nontrivial 
power law type of correlations of the local mass density are found that may extend an order 
of magnitude beyond the fiber length. This implies that there are nontrivial interactions 
present during the formation process. A microscopic understanding of this is currently 
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missing. Because of this, phenomenological deposition models may be useful in studying 
how various effective interactions influence the mass density distribution of the deposit [|12| . 

In addition to their practical applications, random deposition models have been the 
topic of intense study in their own right. In particular, they have been extensively studied 
in the context of continuum percolation theory [P~3"|-fZS|] . These models have included both 
uniformly random networks of various objects as well as some that include hard and soft-core 
interactions between the constituent particles. The quantity of central importance in these 
studies is the percolation threshold or critical particle density which for permeable objects 
can be related to the excluded volume of the particles fl9j . This quantity depends on the 
geometrical shape of the deposited particles as well as on interactions between them in a 
nontrivial way. 

Motivated by processes involving deposition, sedimentation and the theory of continuum 
percolation, we introduce and present a detailed study of the properties of a simple 2D 
deposition model in this work. In the model, the tendency of the deposited particles to 
form clusters is taken into account through a parameter < p < 1. For p — 1, the 
deposition process is uniformly random while for p = 0, only a single connected cluster 
grows from an initial seed fiber |]3U|] . We study the growth dynamics, percolation thresholds 
and spatial correlations of the structures formed by the model in detail, using both analytic 
and numerical methods. In particular, we discuss our results in view of the experimentally 



observed mass density correlations in paper sheets [12 



The outline of this paper is as follows. We first introduce our deposition model in Sec. 
II. In Sec. Ill, we examine the growth properties of single clusters in the limit p = 0. After 
this we present both numerical and analytic results for the percolation threshold of fibers 
and various other objects as a function of p. Section V deals with the theory of spatial 
correlations in networks formed by the model; in particular we concentrate on the form of 
the pair distribution and two point mass density correlation functions. Finally, section VI 
contains our summary and conclusions. 
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II. THE MODEL 



The model studied in this paper, called the "Flocculation Model" is in part motivated 
by the tendency of fibers to adhesively stick during paper-making [0,27]. The model is a 



continuum deposition model defined on a 2D surface. Spatially extended objects, such as 
disks and fibers with a finite area, or widthless needles are individually deposited by the 
following rules. First, both the orientation of the object as well as its spatial coordinates 
are chosen from a uniformly random distribution. If a deposited object lands on another 
object already on the surface, the attempt is always accepted. However, if it lands on empty 
space the attempt is accepted only with a given probability p. Thus, the parameters that 
characterize the model are the acceptance probability p, the dimensions of the deposited 
objects, the linear dimension L of the surface, and the number of deposited objects kept iV 
(i.e. the number of accepted attempts). In this work the deposited objects are mostly fibers 
(rectangles of finite area and varying dimensions); however, in Sec. IV we also consider 
needles (widthless sticks) as well as disks. 

In the limit p = 1 the model reduces to the extensively studied case of a uniformly 
random network p.p^|-|T^, ^ , p^] . However, for p < 1 it is clear that there are effective 
interactions between the particles that tend to enhance cluster formation within the deposit. 
In particular, for the extreme case of p = 0, only particles that always touch each other are 
accepted, and a single, connected cluster grows (assuming there is an initial seed particle). 
We expect all the properties of the growing network, including its percolation threshold 
and spatial mass density correlations to depend on p in a nontrivial way. We note that 
the definition of the model here is different from the extensively studied class of Random 
Sequential Adsorption models, where particle overlap is not allowed [|TJ . 

III. GROWTH OF INDIVIDUAL FIBER CLUSTERS 



First we examine the deposition of fibers in the limit p = of the model |30| where only 



one, connected fiber cluster grows from an initially deposited seed fiber. Fig. 1 shows a 
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typical cluster, for N = 2000 fibers of length A = 1 and width uo = 1/4. An individual 
cluster has an irregular, ragged shape for which an average radius Ri(N) can be defined. 
When a statistical average is taken over all possible clusters of size N, the average cluster 
is spherically symmetric and one can define an average radius R(N) = (Ri(N)) for which 
an expression can be obtained by the following argument. Consider the more general d 
dimensional case of a cluster of iV >> 1 identical fibers of linear dimension a where additional 
AiV fibers are deposited. Its radius grows by AR and area by AA The change in the cluster 
area AA oc a d (R d ~ l / R d )AN oc R d ~ 1 AR 1 which in the continuum limit gives 

R(N) = BN 1/{d+1 \ (1) 

where B is a constant proportional to a. In particular, for the two-dimensional case this 
means that the average size of the cluster grows as iV 1 / 3 as already shown in Refs. [29,130 . 



The radial mass density of the average fiber cluster was also examined in Refs. p9| , p0 
For iV >> 1 (typically N > 500) required for the scaling relation Eq. ([[]), it was shown 
numerically in Ref. JHJ that all the density profiles scale onto one curve. The approximate 



form of this profile can be obtained by the following argument following Ref. Let 



p(r, N) denote the local mass density in a fiber cluster of N fibers, and K the constant flux 
of mass per unit area onto the cluster of radius R(N). Assuming that the rate of change of 
p is proportional to the flux and inversely proportional to the area of the cluster, within the 
cluster radius, then we can write the equation for p as 

For r > R(N), we assume the flux to be zero. Using Eqns. (P and (0) we obtain 

p(r , N) = — R{m - w)) . (3 ) 



Eq. (H) is in good agreement with the numerically obtained density profiles in 30 



However, as already noted in Ref. |J0|, at the edge of the average cluster there is a 



systematic deviation from the linear decrease of the density profile to a slower decay law. 



This indicates that the edge of the cluster is rough, and that the roughness also grows as a 
function of N. This can be characterized by calculating the average width from the definition 

w(N) = [((R t (N)-R(N)) 2 }} 1/2 . (4) 

Since the nonlinear edge of the density profile scales, too, we expect that w(N)/R(N) = 
const, for N » 1. We have numerically simulated w(N) in 2D and indeed find that it 
grows proportional to N 1 ^ 3 . 

IV. PERCOLATION PROPERTIES OF THE MODEL 

When enough fibers are deposited in a finite system, its edges become connected and 
percolation takes place in the model. The properties of the corresponding continuum perco- 



lation transition are of particular interest [19,31 1. For the case of single cluster growth, the 



growth law of Eq. (1) implies that N oc R 3 . Thus, within a finite system L a single cluster 
spans the system only when R « L which means that the density of fibers n = N/L 2 oc L, 
i.e. it increases with the system size. This means that the critical density i] c for percolation 
becomes infinite for p = 0. However, for p > there is a finite probability of nucleation of 
multiple clusters and thus n c should be finite. Indeed, continuum percolation of 2D rectangles 
disks |§05§g§ggigTlgggg and needles ^,0^T7|Q^S| among other geometric 



objects has been extensively studied, and the corresponding critical densities determined 
numerically. In this section, we present results of numerical and analytic calculations of the 
percolation properties of our model for < p < 1 and compare the limit p = 1 with the 
existing studies. 

A. Numerical results 

We first present results of extensive numerical simulations of the critical densities of 
percolation rj c (p) as a function of p for three different geometric objects: needles of length 
A, fibers (rectangles) of length A and width uj (A > u), and disks of radius r^. The details of 
the method we have used are as follows. First, each system contains an inner box of length 
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L and an outer box L + L', with free boundary conditions. It is important in the present case 
to choose V to be large enough so that the average density across the system is constant 
within the inner box f3~2]| . The centers of the objects are distributed within the outer box. 
However, only objects partially or completely within the inner box are allowed to belong to 
the connected clusters. To keep track of the clusters we employ a simple labeling technique 



In analyzing the data we distinguish two different types of percolation processes. In the 
one-sided case, any two opposite edges of the inner box percolate while in the two-sided 
case, all four edges must become connected. In our final analysis, only one-sided percolation 



data was used [|34|| - We obtained the critical density r/ c by the following procedure. We 
deposited a fixed number of objects N for a given system size L and checked for percolation. 
By repeating this 100-1000 times the spanning probability, i.e. the fraction of percolating 
ensembles were obtained. Repeating this procedure for various values of N, the whole curve 
of spanning probabilities was obtained for each system size L, as a function of the density 
1] = N/L 2 . The point where two such curves for any different system sizes intersect, gives an 
approximate value for r\ c . We obtained these points by fitting the curves to error functions. 
The estimates thus obtained were extrapolated and the final values of r] c obtained using 



the standard Monte Carlo renormalization group method [^,^T|], with the smallest system 
studied being the reference system. For completeness, we have also evaluated the correlation 
length exponent v for our model from the MCRG procedure. We note that it is expected to 



have the well-known value of 4/3 as in lattice percolation |19| , pT 



1. Needles 

Typical configurations generated by the model are shown in Figs. 2 (a)-(c) for N = 2000 
needles for various values of p (1.0,0.8,0.6,0.4,0.2,0.05,0.01,0.005, and 0.001). The existence 
of well-defined clusters in the network for the smallest value of p = 0.001 is clearly visible. 
Employing the MCRG procedure for 1000 ensembles and for system sizes L = 6, 20, 30, 40, 60 
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for 0.2 < p < 1.0 and L = 6, 20, 40, 60 for 0.005 < p < 0.05, we obtain the values for rj c {p) 
as displayed in Fig. 3. The curve displays interesting behavior in the two limits p — > 
and p — *■ 1. The expected divergence of r/ c (p) in the limit p —>■ is clearly visible, although 
for smallest p our results here and below are only qualitative due to strongly increasing 
finite size effects. Moreover, r] c (p) approaches the limit p — > 1 approximately linearly with a 
minimum at p ~ 0.2. These results together with the minimum in the curve visible around 
p « 0.2 will be further discussed in Sec. IV B. 

Our best estimate for T] c (l) = 5.64 ±0.02 agrees rather well with other numerical studies 



reported in the literature for the uniformly random case (see, e.g., Ref. ||23||). It also agrees 
with the estimate of Ref. |16| which gives a value of 5.61. Figure 4 shows the effective 
correlation length exponent v(p) as evaluated from the MCRG procedure. We note that 
v(p) is consistent with v = 4/3 in the limit p — > 1. However, for p — > the effective v 
increases. In Sec. V we show that this is due to the increasing range of spatial correlations 
in the case where individual clusters become large. This induces strong crossover effects 



towards long-range correlated percolation type of behavior p6|. 



2. Fibers 

In the case of fibers, we have used objects with an aspect ratio of A/a; = 4. 
Typical networks generated are shown in Figs. 5 (a)-(c) for various values of p 
(1.0,0.9,0.8,0.7,0.5,0.4,0.3,0.2,0.1,0.05,0.01,0.005, and 0.001). Following the procedure em- 
ployed for needles for system sizes L = 20,40,80, we obtain the n c (p) curve as shown in 
Fig. 6. These results have been obtained by averaging over a 100 ensembles. Again, we 
see the expected divergence of r] c (p) as p goes to zero. However, in the other limit p — > 1, 
r] c (p) is almost constant within the error bars. Our best estimate for the p = 1 limit is 
?7 C (1) = 2.73 ± 0.2. To our knowledge, there are no other numerical studies of t] c (l) for the 
present aspect ratio. However, r/ c (l) can be approximately determined by using the excluded 



volume arguments of Ref. [16|, where r] c (l)(A) = const, w 3.57, and the excluded volume 
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(A) = 1.379 [[37]]. This gives results within about 10 % agreement with our numerical data. 
Finally, the correlation length exponent v(p) is shown in Fig. 7. Again, as p — *■ 1, v{jp) is 
consistent with v = 4/3. Moreover, the effective v(p) again increases as p — > 0. 

For completeness, we have also studied the case of isotropic disks. Typical networks of 
N = 1000 disks of radius r d = 0.5 for various values of p (1.0,0.9,0.8,0.7,0.6,0.5,0.45,0.4, 
0.35,0.3,0.25,0.2,0.15,0.1,0.05,0.02,0.01,0.005, and 0.001) are shown in Figs. 8 (a)-(c). Fol- 
lowing the procedure employed for needles for system sizes L = 4, 24, 32, 44, 64 and averaging 
over 800 ensembles for disks of radius = 1, we obtain the n c (p) curve as shown in Fig. 9 
f38fl . Again, rj c (p) diverges as p — > 0. However, for about p > 0.3, r) c (p) ~ const, within the 
error bars. We note that our best estimate for rj c (l) = 0.374 ±0.016 agrees within the errors 
with other numerical studies reported in the literature (see, e.g., Ref. [ZT|). The correlation 



length exponent u(p) is shown in Fig. 10. Again, as p — > 1, v(p) is consistent with v = 4/3. 
Also, the effective v(p) increases as p — ► 0. 

B. Mean-field theory for percolation thresholds 

We can qualitatively understand the behavior of the percolation thresholds for the two 
limits where p — > and p — > 1 using mean- field theory. First, let us discuss the case where 
p — >• 0. To this end, let us define q as the probability of a given object intersecting any other 
object in a uniformly random network (p = 1). This quantity depends on the dimensions of 
the object, and the system size L. For example, in the case of fibers |39 



2(X + uj) 2 2Xu 

q = —^ + ix- (5) 

For needles, one can just set uj = 0. Moreover, for disks a simple calculation gives q = 
Ovrr^/L 2 . 

Next, consider the actual limit of p — > 0. Suppose we begin with an initially empty 
system and drop one object, and then drop another. The probability that it sticks onto 
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the 2D plane but not on the first one is = (1 — q)p. The probability that it does 
stick on top of the first one is = q. Continuing this process for n particles we find 
correspondingly that for the nth particle, = (1 — q) n ~ 1 p, while = 1 — (1 — q) n ~ l . 
Note that for sufficiently small q, we assume that the initially objects tend to land in a 
uniformly random manner. For n greater than some n = n c there comes a point where every 
additional object will encounter as much occupied as empty space. From this point onwards, 
additional objects will essentially tend to build up clusters around the initial n c seed objects. 
By this reasoning, the mean number of clusters n c can be approximately obtained from the 
solution of 

(1 - ff )«e-l p = !_(!_ g)«.-l ( 6 ) 

which gives 

n. = 1 - (7) 
log(l - q) 

We note that this equation is expected to underestimate the true n c due to the fact that due 
to clustering for p < 1, the effective excluded area of the objects will be smaller than that 
given by the bare q in Eq. (|5|). Let us assume that the growth of any one of the n c clusters 
follows Eq. (1) 

R c = BN l J\ (8) 

where iV c is the number of particles in each cluster. Note that 

N = n c N c , (9) 

where iV is the total number of objects in the network. On the average, the system must 
percolate when the clusters begin to overlap, i.e., when 

/J-2fl, (10) 

where p is the density of clusters, p = n c /L 2 . 
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Substituting Eq. (0) into the expression for p and using Eq. (|H]) we obtain, 

-1/2 



(ii) 



Using Eq. (|9[), we obtain 



L I log(l+p) V 1/2 

Equation ([12]) gives us the mean-field percolation threshold defined as r\ c = lim^oo r\ c ^. In 
particular, for the case of fibers we obtain 



Vc 



V 32^° los / (13) 
This result is expected to become more and more accurate in the limit p — > 0, where it 
becomes 



(a + ^) 2 + a^tt n 
* = V w** tp (14) 



The result that r/ c oc yl/p has also been obtained in a somewhat different manner in 
Ref. [p9[| , in accordance with our results. This universal inverse square root divergence onp 
is in qualitative agreement with the numerical data presented in Figs. 3, 6 and 9. However, 
due to the numerical difficulties in determining r] c accurately in this region, we have not 



been able to quantitatively verify Eq. (H). It is interesting to note that in the small p 
limit, where precolation occurs through the growth of well-defined individual clusters, our 
model is somewhat similar to the cooperative lattice filling model where cluster growth leads 



to percolation [filJl . However, in this case there is no divergence of r] c since overlap is not 
allowed. 

Next we consider the other limit of p —>■ 1. Let us again start with a purely uniformly 
random case (p = 1) at its percolation threshold. The probability of any two particles 
intersecting is given by q as before. Take an arbitrary particle from the network of N 
particles. The probability is now given by 

= (1 - qf- 1 . (15) 
11 



Therefore, the number of non-intersecting objects will be, on the average, 

N s = P( N) N. (16) 

The main qualitative idea in the derivation here is the following. The effect of p < 1 is to 
reduce the number of non-intersecting objects due to clustering, but the backbone of the 
percolation cluster itself can be assumed to stay constant for p not too small. This already 
implies that the percolation threshold decreases as p decreases. 

Let us next quantify this argument. Assume that the number single particles removed 
that do not belong to the percolation cluster, is given by 

N r = (l-p)N s . (17) 

Therefore, the number of particles left at the percolation threshold for such p can be written 

as 

N(p) = N(l)-(l-p)N 8 
= N{1)-{1- P )PWN{1), (18) 

where N(l) is the total number of particles within the uniformly random reference system 
with p — 1. For the critical density we thus obtain 

VcAp) = ^(i) - P^VcAV +pP^ n) Vc,lW 

= Ap + B, (19) 

Therefore, we can conclude that t) Cj l(p) increases linearly with p. This is also visible in our 
numerical data for the critical densities, at least for needles and fibers. To examine the 
theory quantitatively, we can show that for the case of fibers, in the limit L — > oo, 

A = ^(1) exp[-77 c (l)( 2(A + a;)2 + 2Xu)]. (20) 

7T 

We have estimated A from the numerical data for needles, and find that the numerical 



and theoretical values agree within about 50 % |42| . The theoretical slope for disks is about 
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an order of magnitude smaller than for needles, and thus it is not surprising that it is not 
discernible in our data of Fig. 9 within the error bars. Also, for the case of fibers the error 
bars are too large to see the slope. 

The analytic derivations of this section can be used to explain the minimum in the i] c (p) 
vs. p curves. Since in the limit p — > 0, t] c (p) must diverge, whereas for p slightly below 
unity rj c (p) must linearly decrease, it follows that a minimum must exist for some < p < 1. 
Physically, the existence of this minimum can be explained by a competition between two 
mechanisms in the model. First, single cluster growth that dominates for small p is rather 
an inefficient way of forming the percolating cluster. On the other hand, uniform random 
filling of the lattice is inefficient, too, since there are many particles that do not belong to 
the percolating cluster. Thus, there is an "optimum" way of forming the cluster at some 
intermediate value of p, where the critical density is minimal. 

V. PAIR DISTRIBUTION AND THE TWO-POINT MASS DENSITY 

CORRELATION FUNCTION 

In this section we begin by writing down a general expression for the two-point mass 
density correlation function for a random network of objects in which there is some effective 
particle-particle interaction. We then generalize this expression for the case of fiber networks. 
This expression involves the pair distribution function for the centers of the objects. This 
function has previously been derived exactly for the case uniformly random network of 
objects P,^0|. We generalize this derivation for our model, and compare the results with 



numerical simulations. We also examine in detail the form of the two-point mass density 
correlation function for the model. 

A. A general expression for the two-point mass density correlation function 



The two-point density correlation function is generally defined as [(43 



G(r) = lim f-\ 1 -^), (21) 



13 



where / is the 2D Fourier transform operator, p(k) is the Fourier transform of the mass 
density distribution and A = L 2 is the area. This distribution can be written as 

N 

p(f) = ^Mf,4^) (22) 
i=i 

where p(x, Xj, 9j) is equal to a over the surface of particle % which is located at the point 
Xi with orientation Oi, and zero outside of it. Here a is the mass density of each particle. 
Writing p{k) = J A p(x)e~ 2mk ' x , and substituting into Eq. (|2ll) we obtain after some algebra 
the general expression |44 

n(k) = ri(\I(k)\ 2 ) + (|/(fc)|) 2 jim N(N 2 ~ 1] (JJ A e-^<*^m-2i\)**i<P*j) > ( 2 3) 

where n(k) = lim^oo \p(k)\ 2 /(Aa 2 ), and r\ is the number density of the particles. This gives 
the correlation function as 

G(r) = r\o*n(k)). (24) 

In the equations above, (\I(k)\ 2 ) denotes the square of the Fourier transform of a single 
object averaged over all orientations. The second term in Eq. ([23|) contains the center- 
center pair distribution function g(\£i — Xj\) Q that describes the effective interparticle 
interactions in the network. It gives the probability per unit area squared of particle center 
i occupying position Xi and that of another particle center j occupying position Xj. Making 
a change of variables and carrying out one of the integrals in Eq. (|23|) gives 

h(k) = v(\I(k)\ 2 ) + 7] 2 (\I(k)\) 2 lim ( [ 2L J (2nkr)[An(r)]dr\ , (25) 

where Jq(x) is the Bessel function of order zero, and we define AQ(r) = lurgir). Equation 
(p5|) has two components. The first describes the self-correlation of single particles, while 
the second describes the cross-correlations, through the function fl(r). The function Q is 
defined as 

. , /# of pairs of centers in a shell (r, r+dr)\ , . 

illrjar = ( — — - — : : ), (26) 

x total jf- of pairs of centers m system ' 
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where the averaging is over all configurations. The function Q(r) in our final expression Eq. 
( p5j) contains the effects of particle interactions, while I(k) contains the information about 
the geometry of the deposited objects. In the special case of fibers, I{k) has been derived 
exactly in Ref. JEJ giving 

1 f 2ir \2 



/ 1 r 2 ' K \ 
(|/(A;)| 2 ) = (\oj) 2 — / sinc[7rfcAcos(t)]sinc[7r£xusin(t)]dt (27) 
\2n Jo ) 



and 



(|/(A;)|) 2 = (Xlo) 2 — [ ^ sine 2 [vrfcA cos(t)]sinc 2 [7r ku sm(t)]dt, (28) 
2n Jo 

where sinc(x) = sin(x)/x. In the previous studies, the second term in Eq. (|25"D has been 
neglected for uniformly random fiber networks. However, for our model, where effective 
interactions arise for p < 1 the second term gives rise to nontrivial correlations that must 
be taken into account. 

B. The pair distribution function 

An exact expression for Q(r) in the case of uniformly random networks has previously 
been derived by Ghosh P|,f40||. Below, we will generalize this derivation to obtain an approx- 
imate formula for our model that is valid for small values of p. To begin with, let ANcM(xi) 
be the number of centers of masses within the area element AA(xi) around the position 
vector Xj. The number of pairs included in two different area elements is then given by the 
product ANcM(xi)ANcM(xj). Further, the total number of pairs separated by vector x in 
a given configuration is 

£ AN C m(x )ANcm(x + x). (29) 

Dividing this by iV(iV — l)/2 we obtain the probability to find a pair of centers of mass 
separated by the vector x in a given configuration. Multiplying and dividing by area elements 
yields 

«*)* - A ^° )A ^ty a w „ + „>. (so, 
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where we have assumed that the system is isotropic and thus Q depends our e \x\ only. 
The brackets denote configuration averaging. Taking the continuum limit where t\cm = 
limAA^o(AA r c < M /AA) we obtain 

2 r 

n(r)dr = ^ _ { J^d 2 x r]cM{xo)T]cM{xo + x)AA(x + £)). (31) 

Hence, 

n(r)dr = 2 _ - jf d 2 x G CM (r , r)AA(r + r), (32) 

where 

Gcm{t q , r) = (t]cm(xo)vcm(xo + %))■ (33) 

It is easy to see that for a uniformly random set of points with translational invariance 
GcM{ r ) — N(N — l)/(2v4 2 ) = const, (no double counting of pairs), which leads to AQ(r) = 
f2/j(r, L) with fl^ being the previously derived exact pair distribution function for a uniform 
random network P,^0|. For large system size (A — > oo), AQ(r) — > lisr. When this is 
substituted into Eq. (^5|) , a constant term arises in the correlation function since there are 
no particle-particle interactions in this special case. 

In the case of p < 1 where clustering occurs in the model, the pair distribution function 
Q(r) gives rise to nontrivial correlations. In this case Q(r) can still be derived approximately 
using mean field arguments. A straightforward but lengthy calculation gives the result (see 
the Appendix for details) 

n(r)=u 1 (N,p)g(r)+u 2 {N,p)n h (r,L), (34) 

where 

f (30r/42i? 5 )(3r 3 - 6Rr 2 + 4R 3 ) , for < r < R; 
g {r) = \ (30r/A2R 5 )(2R - r) 3 , for R < r < 2R; (35) 

[ , for 2R < r. 

The coefficients are defined as LJi(N,p) = n c N c (N c — 1)/[N(N — 1)], where n c and N c are 

defined in Sec. IV, and u 2 (N,p) = 1 — u>i(N,p). In Eq. ( |3~4] ) Q h (r,L) denotes the pair 

distribution function for a uniformly random system that is given by [fl flCH 
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(4r/L 4 )[vrL 2 /2 - 2rL + r 2 /2] , for < r < L; 

Q h (r,L) = { (Ar/L 4 )[L 2 (arcsin(L/r) - arccos(L/r)) (36) 
+2Ly/r* - L 2 - |(r 2 + 2L 2 )} , for L < r < y/2L . 

As discussed in more detail in the Appendix, Eq. (|3^) is only approximately valid, 
although it does reduce to the exact limit £\(r, L) for p — 1. It has been derived under the 
assumption that the deposited objects form clusters that grow independently following the 
scaling relations derived in Sec. III. This approximation is best justified in the p — > limit. 

In Figs. 11(a) and (b) we illustrate several numerically obtained curves for Q(r) for 
p = 0.001 and 0.01, respectively, in the case of fiber networks. There are two common 
features in the figures. First, for small values of N well below percolation, Q(r) develops a 
sharp peak near the origin. This peak is due to the cluster structure in the growing network. 
Moreover, its position is proportional to the average size of the clusters as can be seen from 
Eq. (|35D . In the opposite limit of N — > oo, Q(r) approaches flh(r, L) as expected, since 
the effect of p becomes gradually less important after the network has percolated. In the 
intermediate regions, a combination of the two effects is seen and a smooth transformation 
occurs between the two regimes. We note that when p increases, the height of the first 
peak decreases and its position shifts to smaller values of r indicating that clusters become 
smaller. 

To quantitatively test the accuracy of our analytic formula, we compared the numerically 
obtained function Q(r) with our analytic result of Eq. (|35|). In Figs. 12 (a) through (c) we 
show the results of comparison for p = 0.001, with varying N. In evaluating the analytic 
formula two quantities, namely the average cluster radius R and the number of clusters 
n c are needed. The latter can be evaluated from Eq. (|7|) for any p. However, since this 
equation underestimates the true number of clusters in the system, e.g. for p = 0.001 it 
gives n c pa 1.2, we have allowed n c to be a fitting parameter. The quantity R, on the other 
hand, is fixed using the scaling relation R = BN^ 3 where iV c is the number of fibers in 
each cluster. This can be evaluated from the identity iV = n c N c (see Sec. IV B). The only 
parameter left is the prefactor B which we fix to be B = 0.8. This value has been extracted 
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from numerical data for single cluster growth. As seen in Figs. 12 (a)-(c), setting n c = 4.5 
gives very good agreement between the theory and the simulations. This agreement also 
lends indirect support to the inverse square root divergence of the percolation threshold on 
p derived in Sec. IV B. 

We should note that the derivation of the analytic formula is expected to be valid at 
most up to the point of percolation and for small values of p, since otherwise the whole 
picture of growing independent clusters breaks down. We have not carried out a systematic 
study of the range of validity of the theory for larger values of p, however. 

C. The two-point mass density correlation function 

Next, we examine the two-point mass density correlation function for the model. In 
particular, we examine the case of fiber networks with an aspect ratio of X/u = 20/1, which 
is rather close to the needle limit. We also discuss the effect of different object geometries 
on the correlation function. It can be obtained from our formalism by substituting the 
numerically evaluated Q(r) into Eq. (B^), and using Eq. (E3). However, for simplicity we 
have calculated the correlation function directly from simulations on a lattice ES[ . Defining 



G{x) = {[m(x') — (m)] [m(x'+x) — (m)]) to be the mass density- density fluctuation correlation 
function, we find that its form is well approximated by 

G(r) ~ r ~ aiN > p \ forO < r < A(N,p), (37) 

where A is an effective cutoff for the power law form. In the case of p = 1, we obtain 
a ~ 1, and A rs A/2. For lower values of p, however, G(r) displays interesting behavior, as 
illustrated in Fig. 13, where we show G(r) on a logarithmic scale for p = 0.001. By examining 
the curves it can be inferred that a(N, p) goes through a minimum as A" increases. Moreover, 
A attains a maximum where a is minimum. 

Figure 14(a) shows the effective a(N,p) vs. A" for the model, for different values of p. 
The behaviour of a can be understood as a competition between individual fiber clusters 
and uniformly random fibers, both of which coexist for < p < 1. The approach of the 



correlation exponent toward a = 1 as iV — > oo is caused because in this limit the effect of 
p gradually becomes unimportant as discussed in Sec. V B. On the other hand, the initial 
decrease of the effective a(N,p) for small values of N (as well as the increase of the range 
of the corresponding power law regime) is due to the growth of essentially independent fiber 
clusters. 

A better understanding of the low N behaviour of a(N,p), as well as the increase in 
the power law regime, can be understood by considering the particular case of p = 0. In 
this limit, only a single cluster emerges (see Fig. 1). We find that G(r) as averaged over 
an ensemble of clusters of size N is also well approximated by a power law form. Figure 
15 shows a(N,p = 0) vs. N for this case. For the case of a large enough N (which is 
about 300-400 for the fibers here), a saturates to a fixed value of about 0.05. Moreover, 
we also find that the range of the power law is proportional to the radius of the cluster, 
with a constant of proportionality of about 0.7. We note that a reaches its saturated value 
approximately at the same value of N where the scaling laws of single cluster growth in Sec. 
Ill start to become valid. 

We have also examined the correlation function for the fiber aspect ratio of X/u = 20/5 
that is away from the needle limit, and for isotropic disks. In the case of such fiber or disk 
networks, we find that for values of N and p comparable to those of the needle-like case, the 
correlation function is not as well approximated by a simple power law behavior. However, 
as expected the power law behavior does eventually emerge for small enough values of p. 
The same conclusion also applies to single cluster correlations where larger values of N are 
needed. Thus, we can conclude that the needle-like anisotropy enhances the emergence of 
power law type of behavior in the system. 

It is instructive to consider a random cluster network composed of a uniformly random 
distribution of fiber clusters, each of N fibers, in the needle-like limit. Following the method 
of Sec. V, it is straightforward to show that the correlation function G(r) for such a network 
is given by Eq. (pop, where (|/| 2 ) now represents the square of the Fourier transform of a 
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single cluster. In this case G(r) is exactly proportional to the correlation function of single 
clusters discussed above, and the same power law form emerges, with its range proportional 
to the cluster size. However, we note that as discussed in Ref. | 5Cfl , even longer range power 



law type of correlations for a network of fiber clusters can be obtained by assuming that 
there is an effective cluster-cluster interaction that gives rise to such correlations beyond the 
cluster size. 

It is interesting to examine our results for the areal mass density fluctuation correlation 
function with recent measurement of this quantity for laboratory made paper fll2| . It has 



been shown that G(r) for these paper sheets is well approximated by a power law form 
whose exponent is independent of the mass density of paper, and is a ~ 0.37. Moreover, the 
range of this power law is rather short for high basis weight paper, i.e. of the order of the 
fiber length, whereas it becomes longer ranged for low mass density paper extending even 
up to about 14 times the fiber length. It is interesting to note that the simple deposition 
model presented here can give similar values of a for values of N corresponding to real 
paper, by suitably adjusting p (see Fig. 14). However, the range of the power law in this 
case does not vary with mass density in the same way as in real paper. We can conclude 
from the results of this work that there are at least two possibilities of obtaining longer 
range correlations such as in low mass density paper. The first is that paper comprises small 
clusters between which there is an effective interaction giving rise to such correlations as 



explained in Ref. |P0(1 . The second is that paper can be described by deposition of larger, 



independent clusters whose internal correlations are of power law type, with a ~ 0.37 |L2 
However, the true microscopic modeling of paper-making and unraveling of the reasons 
behind spatial correlations is a formidable theoretical problem. 

VI. SUMMARY AND CONCLUSIONS 



In this paper we have studied growth, percolation and correlation properties of disordered 
networks of fibers and other extended objects. In particular, we have introduced a deposition 
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model that introduces an effective sticking of particles upon deposition. The key parameter 
of the model is the acceptance probability p which describes the effective interactions between 
particles, thereby controlling the tendency to form clusters. This choice of rule is motivated 
by interactions that may be present in suspension during the formation phase of paper. 

By varying p the model interpolates between the growth of uniformly random networks 
(p = 1) and single cluster growth {p = 0). We started by studying the growth of single 
clusters, and the form of the density profiles for such clusters. For both quantities, analytic 
results have been derived and verified by numerical simulations. Following this, we examined 
the percolation threshold i] c (p) of needles, fibers, and disks in detail by performing extensive 
numerical calculations as a function of p. In the limit p — 1, our results reduce to those 
previously obtained for needles and disks. In the case of rectangles, however, we were 
unable to make quantitative comparisons with the literature. Furthermore, our results were 
supplemented by mean field arguments for the dependence of r] c {p) in the limits p — > 1 and 
p -> 0. 

An additional interesting feature of the model is the emergence of nontrivial spatial 
correlations for p < 1. We have examined these correlations by an approximate analytic 
calculation of the pair distribution function of a network of objects, and have numerically 
examined this function for the case of fibers. Moreover, we have studied the two-point mass 
density fluctuation correlation function of the model, and shown that in some cases it can 
be well approximated by a power law type of form, with a nontrivial exponent. This kind 
of power law behavior is enhanced by increasing anisotropy of the deposited objects. 

Finally, the results for the correlation function from the model are particularly interesting 
when contrasted with those obtained for mass density correlations in real paper sheets, 
where power-law type of correlations arise in some cases. Further work on related, more 
complicated deposition models for fibers is in progress. 
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APPENDIX A: THE PAIR DISTRIBUTION FUNCTION 

In order to derive an expression for Q(r) we need to know how to perform the ensemble 
average in Eq. fl33l). To derive an explicit form for Gcm{xq,x) we start with a general 
expression for the probability density P[p] of a given fiber configuration, where p is the mass 
density: 

. N N 

p ip] = / II d 2 x k B{{x k , 6 k }) 6[p(x) - $s, 9.)) , (Al) 

JA k=l s=l 

where p has been defined in Eq. (ESj) The probability distribution for the centers of masses 
is readily obtained by replacing p by a delta function. 

To calculate correlation functions the distribution function of the centers of masses B 
has to be provided. In what follows, we will for simplicity drop the angular dependence from 
it. In case we assume that p is so small that during deposition single clusters emerge and 
start growing independently of each other, the distribution will factorize as 

B({x k }) = 1] B cl ({x s } m ) , (A2) 

m=l 

where {x s } m denotes the set of the center of mass coordinates of the particles belonging to 
the m th cluster, and n c is the number of the clusters. The assumption made here is that 
each particle must be uniquely assigned to one of the connected clusters in the network. The 
distribution of the particles inside each cluster B c \ is taken to be a completely symmetric 



function of all of its arguments. In Eq. (33) we will separate the different terms in the 



correlation function into two groups. First, there are those terms that come from correlations 
between particles belonging to the same cluster (corresponding particle labels are i and j). 
The remaining contribution is collected from correlations between particles sitting in different 
clusters (labels I and J): 

N 

GcM(x ,x)=y2 IT d 2 x k B({x k }) 8(x - Xi) 8(x + x- £,) (A3) 



i,j JA k=l 
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r N 

+ H / II d 2 x k B({x k }) 8(x - xi) 8(x + x- xj) (A4) 
i j JA k=i 



Due to the indistinguishability of the particles the first sum on the RHS of the previous 
equation can be replaced with a factor uj\ = n c N c (N c — l)/2. Similarily, the second sum 
gives rise to a factor uj 2 = N(N — l)/2 — n c N c (N c — l)/2. Hence, 

Gcm(xo,x) = Ui / Y[d 2 x p B d ({x p }) 5(x - Xi)5(x + x- xj) (A5) 
+ uj 2 / JJd^pBoj^fp}) 5(S -xj) / 5(x + f-fj) (A6) 

p JA q 

= G) X bV> (f o, x + f ) + cu 2 4 2) (f )^ 2) (f o + f ) . ( A7) 

In the equation above, particle coordinates {x p } belong to one of the n c identical clusters 
(xi,Xj e {x p }), and x q belongs to another cluster (xj G and xj G 

The quantity B d ({x p }) \~\ p d^Xp tells the probability to find the centers of the masses of the 
particles within one cluster at positions x p . Relative to the center of mass x c i of the cluster 
these centers of masses cannot be scattered around the whole system at random positions 
for the cluster is assumed to form a connected structure. From single cluster growth theory 
we know that the mass density p is a cone-like function (given by Eq. (H)) with radius R. 
We assume that the number density ?]cm of the centers of masses of particles is proportional 
to the mass density p in the cluster. The first term B^\xq, xq + x) is obtained by integrating 
out the positions of all the other centers of masses of particles but two, which are fixed at 
the positions Xq and Xq + x. A nonzero contribution from the integrations will emerge only 
from those regions where the two fixed points are within a distance R from the center of 
mass x c i of the cluster, i.e. when they belong to the cluster. Furthermore, we assume that 
the number of particles in the cluster is large enough so that x c \ is independent of the two 
points xo and x*o + x. Thus, we can write 

fl«(Sb, % + g) = Cl ! A ci p(f - swrfib + i - Sa) ; (A8) 

J A 

B % {&o) =C 2 d 2 x d p{x - x d ) , (A9) 
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where C\ and C2 are normalization constants (B c i is a probability density). The integration 
domain can be extended over the whole system since p(r) = when r > R. Substitution of 
Gcm{x , x) into Eq. fl32|) yields 



Q(r)dr = uj\g(r)dr + u^f^r, L)dr , (A10) 
where u x = Q 1 ^2 = ^ iwfiT , and 



g ( r )dr= / diro-By(fo,^o + ^)AA(fo + ^) (All) 

JA 

= 2nrdrCiL 2 / d 2 xop(x*o — x c i)p(xo + x — x c i); (A12) 
Vt h {r,L)dr = C\ \ d 2 x / d 2 x c ip(x — x c {) \ d 2 x d p(x + x — x* cl )AA(x + x) (A13) 

JA JA JA 

= / d 2 x AA(fo + J?). (A14) 

JA 

where the integration domain K in Eq. ( A12 ) is given by \xq— x c {\ < R, and \xq+x— x c i\ < R, 



which is the region where there in a nonzero contribution to g(r). This means that we can 
make the transformation of variables u = Xq — x c i, 2v = Xq + x c i, and thus the final result 
only depends on r = \x\. 

In d = 2 d = (LnR 2 a/3)~ 2 and C 2 = (L 2 nR 2 a/3y 1 with a = p(0,N). The density p 
is given in Eq. (^|). Plugging it into Eq. ( |A12j ) and setting AA(xo + x) = 2nrdr, which 
is true within each cluster, yields g(r). To get some insight how g behaves we will use one 
dimensional density profiles instead of two dimensional ones. However, we do not replace 
AA(x + x) with its one dimensional counterpart, AA(xq + x) = 2 dr since we want to satisfy 
the relation g(0) = 0, which is true in d — 2. In this case, the normalization constant C\ is 
given by C\ = 15/(7L 2 irR 3 a 2 ). This gives as the final result the Eq. (|35|). 

Finally, we would like to note that in Eq. (|A14|) AA(xq + x) imposes a cutoff at the 



edges of the system. This area element can be written in the form 

d 2 x AA(x + x) = rdrfd 2 x P(x ,x), (A15) 



where (3(xo,x) tells the maximum angle swept by a vector of length r = |x| rotated around 
a fixed position xo, with the constraint that the vector must lie inside the system. This can 
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be used in deriving the pair distribution function of a uniformly random system as given 
by Eq. pi). 
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FIGURE CAPTIONS 

Fig 1. A snapshot of a single growing fiber cluster for N = 2000 fibers of length A = 1 and 
width uj = 1/4. 

Fig 2. Snapshot of a disordered network of N = 2000 deposited needles of length A = 1, for 
(a) p = 1.0, (b) p = 0.1, and (c) p = 0.001. 

Fig 3. The critical percolation threshold rj c (p) vs. p for a network of needles of length A = 1, 
extrapolated to L — > oo. The error bars are of the size of the points. 

Fig 4. Correlation length exponent v(jp) vs. p for needles of length A = 1. In this and the 
following figures for v(p), the horizontal line denotes the exact value v = 4/3. 

Fig 5. Snapshot of a disordered network of N = 2000 fibers of length A = 1 and width 
uj = 1/4. In (a) p = 1.0, in (b) p = 0.1, and in (c) p = 0.001. 

Fig 6. The critical percolation threshold rj c (p) vs. p for a network of fibers of length A = 1 
and width uj — 1/4, extrapolated to L — > oo. See text for details. 

Fig 7. Correlation length exponent v(p) vs. p for fibers of length A = 1 and width u — 1/4. 

Fig 8. Snapshot of a network of iV = 1000 disks of radius r d = 0.5 in a system of size L = 20. 
In (a) p = 1.0, in (b) p = 0.1, and in (c) p = 0.001. 

Fig 9. The critical percolation threshold r] c {p) vs. p for a network of disks of radius r d = 1, 
extrapolated to L — > oo. 



29 



Fig 10. Correlation length exponent is(p) vs. p for disks of radius r d = 1. 

Fig 11. Simulated pair distribution functions for fibers of length A = 1 and width 
lu = 1/4 for (a) p = 0.001 and (b) p = 0.01. The different curves correspond to 
N = 15,25,75,150,250,500,750,1000,2000,4000,4500,5000,6000 from top to bottom at 
r = 0. 

Fig 12. Comparison between analytic (solid lines) and simulated (dotted lines) pair distri- 
bution functions for fibers of length A = 1 and width uj — 1/4, for p = 0.001 and L = 20. In 
(a) N = 250, (b) N = 500, and (c) N = 750. See text for details in this and the following 
figures. 

Fig 13. A plot of ln(G(r)) vs. ln(r) for fibers of length A = 20 and width u = 1, for p = 0.001. 
From bottom to top, N increases with N = 25, 50, 300, 1500, 3000, 5000, 8000, 20000, 25000. 
The system size L = 400. 

Fig 14. A plot of the effective exponent a(N,p) vs. iV for various p. The aspect ratio is 
L/\/uj = 400/20/1. 

Fig 15. A plot of a(N,p = 0) vs. N for the special case where a single cluster 
grows with fibers of length A = 20 and width uj — 1 (see Fig. 1). Note the satu- 
ration to a constant a 0.05 after iV 500. The inset shows ln(G(r)) vs. ln(r) for 
N = 25, 100, 300, 600, 1500, 2000, 3000, 4000, 5000, from bottom to top. 
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